home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / ali.cc next >
Encoding:
Text File  |  1995-12-20  |  8.0 KB  |  265 lines  |  [TEXT/ALFA]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *             Aitken-Lagrange interpolation
  7.  *
  8.  * This package allows one to interpolate a function value for a given
  9.  * argument value using function values tabulated over either uniform or
  10.  * non-uniform grid. The latter is specified by a vector of node point
  11.  * abscissae. The uniform grid is specified by a grid mesh and the
  12.  * abscissa of the first grid point.
  13.  *
  14.  * Synopsis
  15.  *    double ali(q,x0,s,y)
  16.  *    double q             Argument value specified
  17.  *    double x0             Abscissae for the 1. grid point
  18.  *    double s             Grid mesh, >0
  19.  *    VECTOR y            Vector of function values
  20.  *                    tabulated at points
  21.  *                    x0 + s*(i-y.q_lwb()))
  22.  *                                      The vector must contain at 
  23.  *                    least 2 elements
  24.  *
  25.  *    double ali(q,x,y)
  26.  *    const double q            Argument value specified
  27.  *    const VECTOR x            Vector of grid node abscissae
  28.  *    const VECTOR y            Vector of function values
  29.  *                    tabulated at points x[i]
  30.  *                                      The vector must contain at 
  31.  *                    least 2 elements
  32.  * Output
  33.  *    Both functions return the interpolated function value at point q
  34.  *    Interpolating process finishes either
  35.  *        - if the difference between two successive interpolated
  36.  *          values is absolutely less than EPSILON
  37.  *        - if the absolute value of this difference stops
  38.  *          diminishing
  39.  *        - after (N-1) steps, N being the no. of elements in vector y
  40.  *
  41.  * Algorithm
  42.  *    Aitken scheme of Lagrange interpolation
  43.  *    Algorithm is described in the book
  44.  *        "Mathematical software for computers", Institute of
  45.  *         Mathematics of the Belorussian Academy of Sciences,
  46.  *         Minsk, 1974, issue #4,
  47.  *         p. 146 (description of ALI, DALI subroutines)
  48.  *         p. 180 (description of ATSE, DATSE subroutines)
  49.  *    The book essentially describes IBM's SSP package.
  50.  *
  51.  * $Id$
  52.  *
  53.  ************************************************************************
  54.  */
  55.  
  56.  
  57. #include "LinAlg.h"
  58. #include "math_num.h"
  59. #include <std.h>
  60. #include <float.h>
  61.  
  62.  
  63. //#define DEBUG
  64.  
  65. /*
  66.  *------------------------------------------------------------------------
  67.  *             Class that handles the interpolation
  68.  */
  69.  
  70. class ALInterp
  71. {
  72.   Vector arg;            // [1:n] Arranged table of arguments
  73.   Vector val;            // [1:n] Arranged table of function values
  74.   double q;                   // Argument value the function is to be
  75.                 // interpolated at
  76.  
  77. public:
  78.                 // Construct the arranged tables for the
  79.                 // uniform grid
  80.   ALInterp(const double q, const double x0, const double s, const Vector& y);
  81.                 // Construct the arranged tables for the
  82.                 // non-uniform grid
  83.   ALInterp(const double q, const Vector& x, const Vector& y);
  84.  
  85.   double interpolate(void);    // Perform actual interpolation
  86. };
  87.  
  88. /*
  89.  *------------------------------------------------------------------------
  90.  *
  91.  *        Arranging data for the Aitken-Lagrange interpolation
  92.  *
  93.  * Abscissae (arg) and ordinates (val) of the grid points should be arranged
  94.  * in such a way that the distance abs(q-arg[i]) would increase as i
  95.  * increases. 
  96.  * Here q is the point the function is to be interpolated at.
  97.  *
  98.  */
  99.  
  100.                 // Construct the arranged tables for the
  101.                 // uniform grid
  102. ALInterp::ALInterp(const double _q, const double x0,
  103.            const double s, const Vector& y)
  104.     : arg(y.q_no_elems()), val(y.q_no_elems()), q(_q)
  105. {
  106.   const int n = y.q_no_elems();
  107.   assure( n > 1, "Vector y (function values) must have at least 2 points");
  108.   assure( s > 0, "The grid mesh has to be positive");
  109.  
  110.             // First find the index of the grid node which
  111.             // is closest to q. Assign index 1 to this
  112.             // node. Then look at neighboring grid nodes
  113.             // and assign indices to them
  114.             // (kind of breadth-first search)
  115.   int js = (int)( (q-x0)/s + 1.5 );    // Index j for the point x0+s*j
  116.                     // which is the closest to q
  117.   if( js < 1 )                // Check for the case of extrapolation
  118.     js = 1;                // to the left end
  119.   else if( js > n )
  120.     js = n;                // or to the right end
  121.  
  122.   register int dir;            // Direction to the next closest
  123.                     // to q grid node
  124.   dir = ( q > x0 + (js-1)*s ? 1 : -1 );
  125.  
  126.   register int jcurr = js, jleft = js, jright = js;
  127.   register int i;
  128.   for(i=arg.q_lwb(); i<=arg.q_upb(); ++i) // Pick up elements x0+s*i
  129.   {                    // in the neighborhood of q
  130.     arg(i) = x0 + (jcurr-1)*s;
  131.     val(i) = y(jcurr-1+y.q_lwb());    // Once the closest to q point js
  132.     if( jright >= n )            // is found, we pick up points
  133.       dir = -1;                // alternatively to the right
  134.     if( jleft <= 1 )            // and to the left of the js
  135.       dir = 1;                // further and further
  136.     if( dir > 0 )
  137.     {
  138.       jcurr = ++jright;
  139.       dir = -1;
  140.     }
  141.     else
  142.     {
  143.       jcurr = --jleft;
  144.       dir = 1;
  145.     }
  146.   }
  147. }
  148.  
  149.  
  150. static inline int fsign(const float f)        // Return the sign of f
  151. { return f < 0 ? -1 : f==0 ? 0 : 1; }
  152.  
  153.                 // Construct the arranged tables for a
  154.                 // non-uniform grid
  155. ALInterp::ALInterp(const double _q, const Vector& x, const Vector& y)
  156.     : arg(x.q_no_elems()), val(y.q_no_elems()), q(_q)
  157. {
  158.   assure( y.q_no_elems() > 1,
  159.         "Vector y (function values) must have at least 2 points");
  160.   are_compatible(x,y);
  161.  
  162.                     // Selection is done by sorting x,y arrays
  163.             // in the way mentioned above. Fisrt an array
  164.                     // of indices is created and sorted, then arg,
  165.             // val arrays are filled in using the sorted indices
  166.   class index_permutation
  167.   {
  168.     struct El { int x_ind; float x_to_q; };    // x_to_q = |x[x_ind]-q|
  169.     El * permutation;
  170.     const int n;
  171.     static int comparison_func(const void * ip, const void * jp)
  172.       { return fsign(((const El*)ip)->x_to_q - ((const El*)jp)->x_to_q); }
  173.   public:
  174.     index_permutation(const double q, const Vector& x) :
  175.     permutation(new El[x.q_no_elems()]), n(x.q_no_elems())
  176.     {
  177.       register El * pp = permutation;
  178.       for(register int i=x.q_lwb(); i<=x.q_upb(); i++,pp++)
  179.         pp->x_ind = i, pp->x_to_q = fabs(q-x(i));
  180.                     // Sort indices so that
  181.                 // |q-x[x_ind[i]]| < |q-x[x_ind[j]]|
  182.                 // for all i<j
  183.       qsort(permutation,n,sizeof(permutation[0]),comparison_func);
  184.     }
  185.     ~index_permutation(void) { delete permutation; }
  186.                     // Apply the permutation to x to get arg
  187.                     // and to y to get val
  188.     void apply(Vector& arg, Vector& val, const Vector& x, const Vector& y)
  189.     {
  190.       register const El* pp = permutation;
  191.       for(register int i=arg.q_lwb(); i<=arg.q_upb(); i++,pp++)
  192.         arg(i) = x(pp->x_ind), val(i) = y(pp->x_ind);
  193.       assert(pp==permutation+n);
  194.     }
  195.   };
  196.   
  197.   index_permutation(q,x).apply(arg,val,x,y);
  198.  }
  199.  
  200. /*
  201.  *------------------------------------------------------------------------
  202.  *            Aitken - Lagrange process
  203.  *
  204.  *  arg and val tables are assumed to be arranged in the proper way
  205.  *
  206.  */
  207.  
  208. double ALInterp::interpolate()
  209. {
  210.   register double valp = val(1);    // The best approximation found so far
  211.   register double diffp = DBL_MAX;    // abs(valp - prev. to valp)
  212.   register int i,j;
  213.  
  214. #ifdef DEBUG
  215.   arg.print("arg - interpolation nodes");
  216.   val.print("Arranged table of function values");
  217. #endif
  218.             // Compute the j-th row of the Aitken scheme and
  219.             // place it in the 'val' array
  220.   for(j=2; j<=val.q_upb(); j++)
  221.   {
  222.     register double argj = arg(j);
  223.     register REAL&  valj = val(j);
  224.     for(i=1; i<=j-1; i++)
  225.       valj = ( val(i)*(q-argj) - valj*(q-arg(i)) ) / (arg(i) - argj);
  226.  
  227. #ifdef DEBUG
  228.     message("\nval(j) = %g, valp = %g, arg(j) = %g",valj,valp,argj);
  229. #endif
  230.  
  231.     register double diff = fabs( valj - valp );
  232.  
  233.     if( j>2 && diff == 0 )          // An exact result has been achieved
  234.       break;
  235.  
  236.     if( j>4 && diff > diffp )        // Difference stoped diminishing
  237.       break;                // after the 4. step
  238.  
  239.     valp = valj;  diffp = diff;
  240.   }
  241.  
  242.   return valp;
  243. }
  244.  
  245.  
  246.  
  247. /* 
  248.  *=======================================================================
  249.  *                Root modules
  250.  */
  251.  
  252.                 // Uniform mesh x[i] = x0 + s*(i-y.lwb)
  253. double ali(const double q, const double x0, const double s, const Vector& y)
  254. {
  255.   ALInterp al(q,x0,s,y);
  256.   return al.interpolate();
  257. }
  258.  
  259.                 // Nonuniform grid with nodes in x[i]
  260. double ali(const double q, const Vector& x, const Vector& y)
  261. {
  262.   ALInterp al(q,x,y);
  263.   return al.interpolate();
  264. }
  265.